Warping and Precession of Accretion Disks Around Magnetic 

Stars: Nonlinear Evolution 



Harald P. Pfeiffer^ and Dong Lai 

Center for Radiophysics and Space Research, Cornell University, Rhaca, NY 14853 
^ . Email: harald, dong@astro.cornell.edu 

O 

o 

(N : ABSTRACT 

The inner region of the accretion disk around a magnetized star (T Tauri 
^ ■ star, white dwarf or neutron star) is subjected to magnetic torques that induce 

warping and precession of the disk. These torques arise from the interaction 
^ ■ between the stellar field and the induced electric currents in the disk. We carry 

^ ■ out numerical simulations of the nonlinear evolution of warped, viscous accretion 

■ disks driven by the magnetic torques. We show that the disk can develop into a 



o 



highly warped steady state in which the disk attains a fixed (warped) shape and 



. precesses rigidly. The warp is most pronounced at the disk inner radius (near the 

. magnetosphere boundary). As the system parameters (such as accretion rate) 

p • change, the disk can switch between a completely flat state (warping stable) and 

Q . a highly warped state. The precession of warped disks may be responsible for 

^ . a variety of quasi-periodic oscillations or radiation flux variabilities observed in 

. many different systems, including young stellar objects and X-ray binaries. 
> ' 

^ . Subject headings: accretion, accretion disks - stars: neutron - stars: magnetic 

. fields - stars: pre-main-sequence - binaries: general 



1. Introduction 



Warped accretion disks are believed to exist in a variety of astrophysical systems, in- 
cluding X-ray binaries, young stellar objects (YSOs) and active galactic nuclei. For example, 
the well-known 164 day precession of the jet in SS 433 and the long-term variabilities ob- 
served in many X-ray binaries (e.g. Her X-1, LMC X-4; see Priedhorsky Sz Holt 1987; Scott 
et al. 2000; Ogilvie & Dubus 2001) have been explained by the precession of tilted accretion 
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disks. The changes in the flow directions of several YSO jets (e.g. Bate et al. 2000) and 
the photometric variabihties of some T Tauri stars (e.g., AA Tauri; see Bouvier et aL 1999; 
Carpenter et aL 2001) have also been associated with warped, precessing disks. 

Several mechanisms for exciting/maintaining warps in accretion disks have been pro- 
posed in recent years. Schandl & Meyer (1994) considered a warping instabihty due to 
irradiation-driven disk wind. Pringle (1996) showed that even without wind loss, radiation 
pressure itself can induce warping in the outer region of the disk (see also Maloney, Begelman 
& Nowak 1998; Wijers & Pringle 1999; Ogilvie & Dubus 2001). QuiUen (2001) showed that 
a wind passing over the disk surface may induce warping via Kelvin-Helmholtz instability. 
Finally, in the case of disk accretion onto magnetic stars (e.g., neutron stars, white dwarfs 
and T Tauri stars), the stellar magnetic field can induce disk warping and precession (Lai 
1999; see also Terquem & Papaloizou 2000). In this paper we are concerned with these 
magnetic effects on accretion disks. 

The study of disk accretion onto magnetic stars has a long history (e.g., Pringle & Rees 
1972; Ghosh & Lamb 1979; Anzer & Borner 1980, 1983; Lipunov & Shakura 1980; Wang 
1987, 1995; Aly & Kuijpers 1990; Spruit & Taam 1993; Shu et al. 1994; van Ballegooijen 
1994; Lovelace et al. 1995,1999; Campbell 1997; Lai 1998; see also numerical simulations 
by Hayashi et al. 1996; Miller & Stone 1997; Goodson et al. 1997; Fendt k Elstner 2000; 
Romanova et al. 2003). Most previous studies have, for simplicity, adopted the idealized 
geometry in which the magnetic axis, the spin axis and the disk angular momentum are 
aligned. However, it was shown in Lai (1999) that under quite general conditions, the 
stellar magnetic field can induce warping in the inner disk and make the disk precess around 
the spin axis (see §2; see also Aly 1980; Lipunov & Shakura 1980; Terquem & Papaloizou 
2000). Such magnetically driven warping and precession open up new possibilities for the 
dynamical behaviors of disk accretion onto magnetic stars. Shirakawa & Lai (2002a) studied 
these effects in weakly magnetized accreting neutron stars and showed that the magnetic 
warping/precession effects may explain several observed features of low-frequency quasi- 
periodic oscillations in low-mass X-ray binaries. Shirakawa & Lai (2002b) also studied the 
linear, global warping/precession modes of inner disks of highly magnetized {B ~ 10^^ G) 
neutron stars (NSs), as in accreting X-ray pulsars, and suggested that magnetically driven 
disk warping and precession can explain the mHz variabilities observed in X-ray pulsars. 

The studies by Shirakawa & Lai (2002a,b) were restricted to linear warping. Such linear 
analysis is only the first step toward understanding the observational manifestations of the 
magnetic warping/ precession effects. However, as we show in this paper, under many circum- 
stances, the disk warping becomes nonlinear, and the dynamics of nonlinear warped disks 
is much richer. This is the subject of this paper. After briefly reviewing the basic physics 
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of magnetically driven disk warping and precession (§2), we present the dynamical equa- 
tions for the warped disk (§3) and our numerical results (§4), and discuss their astrophysical 
implications (§5). The appendix contains more details of our numerical method. 



2. Magnetically Driven Disk Warping/Precession 

For completeness, we briefly review the basic physics of magnetically driven warp- 
ing/precession and give the formulae to be used in later sections. 

The inner region of the accretion disk onto a rotating magnetized central star is subjected 
to magnetic torques that induce warping and precession of the disk (Lai 1999). These 
magnetic torques result from the interactions between the accretion disk and the stellar 
magnetic field. Depending on how the disk responds to the stellar field, two different kinds 

of torque arise: (i) If the vertical stellar magnetic field penetrates the disk, it gets twisted 
by the disk rotation to produce an azimuthal field AB^ = ^(B^ that has different signs 
above and below the disk {( is the azimuthal pitch of the field line and depends on the 
dissipation in the disk), and a radial surface current Kf results. The interaction between 
Kr and the stellar i?^ gives rise to a vertical force. While the mean force (averaging over 
the azimuthal direction) is zero, the uneven distribution of the force induces a net warping 
torque which tends to misalign the angular momentum of the disk with the stellar spin axis, 
(ii) If the disk does not allow the vertical stellar field (e.g., the rapidly varying component of 
Bz due to stellar rotation) to penetrate, an azimuthal screening current will be induced 
on the disk. This interacts with the radial magnetic field Bj. and produces a vertical 
force. The resulting precessional torque tends to drive the disk into retrograde precession 
around the stellar spin axis. 

In general, both the magnetic warping torque and the precessional torque are present. 
For small disk tilt angle /? (the angle between the disk normal and the spin axis), the 
precession angular frequency and warping rate at radius r are given by (Lai 1999) 

where /j, is the stellar magnetic dipole moment, 6 is the angle between the magnetic dipole 
axis and the spin axis, f2(r) is the orbital angular frequency, and S(r) is the surface density 
of the disk. The dimensionless function D(r) is given by 



D{r) = max (^^r^rl-l, ^2H{r)/ri^ 



(3) 
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Here H{r) is the half-thickness of the disk, and r^n is the inner disk radius, given by the 
magnetosphere radius 

'■•""''(^) • 

where M is the stellar mass, M is the mass accretion rate, and r) ~ 0.5 — 1. The function 
F{9) depends on the dielectric property of the disk. We can write 

F{e) ^2/008^9-8111^6, (5) 

so that F{9) = — sin^ 9 if only the spin-variable vertical field is screened out by the disk 
(/ = 0), and F{9) = 3 cos^ — 1 if all vertical field is screened out (/ = !). In reality, / lies 
between and 1. For concreteness, we shall set F{9) = — sin^6' in the following. 



3. Dynamical Equations for Warped Disks 

In a viscous accretion disk (with the dimensionless viscosity parameter a > H/ r) , the 

dynamics of the disk warp is dominated by viscous stress (rather than bending waves; see 
Papaloizou & Lin 1995; Terquem 1998), and can be studied using the formalism of Papaloizou 
& Pringle (1983) (see also Pringle 1992; Ogilvie 1999; Ogilvie & Dubus 2001). We model the 
disk as a collection of rings which interact with each other via viscous stresses. Each ring 
at radius r has the unit normal vector l(r, t). In the Cartesian coordinates, with the z-axis 
along the stellar spin a), we can write 



^sin j3 cos 7^ 
sin (3 sin 7 
cos 13 



(6) 



with P{r,t) the tilt angle and ^{r,t) the twist angle. The relevant equations are mass 
conservation. 
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Here ui is the usual disk velocity (measuring the r-0 stress), z/2 denotes the viscosity in the 
r-z stress which is associated with reducing the disk tilt and N represents external torques 
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acting on the disk. The radial velocity was derived by Pringle (1992), 

rS|: (r2Q) 2 ^ir^) 

The nonlinear warp equation (8) is based on a phenomenological description of the viscosities 
z/i and 1^2 and is formally valid only in the linear regime when the local disk warp (91/(9 In r is 
small (Ogilvie 1999). A more rigorous treatment by Ogilvie (1999) removes this assumption 
and finds that ui and 1/2 effectively depend on the warp amplitude dl/dlnr. Ogilvie (1999) 
also finds an additional term causing dispersive propagation and precession. However, it 
turns out that this latter term is not very important for Keplerian disks. For the current 
investigation into the qualitative features of the magnetic warping and precession instability 
we adopt the simpler equations given by Pringle (1992). 

The torque N in Eq. (8) was derived in Lai (1999) as 
N 

^ ^ilpCosP (2j xl-Tw cos P[(2; - {u ■ 1)1] (10) 

where 0,p and F^j, are given by Eqs. (1) and (2). The first term of Eq. (10) causes retrograde 
precession of the tilted disk (since Qp < due to the choice F{9) — — sin^^), whereas the 
second term tries to warp the disk. 

The magnetic torque formulae [eqs. (1) and (2)] contain uncertain parameters (e.g., C, 
which parametrizes the amount of azimuthal twist of the magnetic field threading the disk) ; 
this is inevitable given the complicated nature of magnetic field - disk interaction (see Lai 
1999 and references therein). Also, while the expression for the warping torque [eq. (2)] 
is formally valid for large disk warps, the expression for the precession torque was derived 
under the assumption that the disk is locally flat [eq. (1) is strictly valid only for a completely 
flat disk; see Aly 1980]; when this assumption breaks down (i.e., when dl/dlnr is large), 
we expect a similar torque expression to hold, but with modified numerical factors (e.g. the 
function D{r) in eq. (1) will be different). When an almost fiat disk becomes unstable to 
the warping instability, and also in the deeply nonlinear warp regime with disk tilted to 
large radn (see §4), the condition \dl/dlnr\ < 1 is well satisfied. Thus we believe that our 
simplifying assumptions capture the qualitative behavior of accretion disks subject to the 
magnetic torques. 



m 

dr 



(9) 



3.1. Boundary Conditions 



We need boundary conditions for 1 and for the surface density E. The torque (10) decays 
rapidly with radius, so that we expect the outer disk to be undisturbed. This is confirmed 
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by the linear stability analysis of Shirakawa & Lai (2002a,b). We therefore adopt simple 
Dirichlet boundary conditions at the outer edge of the disk, 

krout) = ^out = const. (11) 

Throughout this paper we choose lout = <-^, such that the outer disk is perpendicular to the 
rotation axis of the star. At the inner edge of the disk, we use a zero torque boundary 
condition, 

|(nn) = 0. (12) 

The boundary condition on S is motivated by the flat Keplerian disk, which has^ 

^flat = ^ J , (13) 

with J — 1 — (1 — Jin) (r/rin)^^^^. We choose Dirichlet boundary conditions on the surface 
density, 

M M 

S(^out) — 7^ >^outj S(^m) — 7^ >Jim (14) 

where Jout = 1 — (1 — Jin){T out / fin)"^^'^ ■ We take Jin to be free, parametrizing the physics 
at the inner edge of the disk. 



3.2. Dimensionless Equations 

The radius is measured in multiples of the radius of the inner edge, x — r/rim and we 
introduce a dimensionless surface density cr = E/Eq. We assume viscosities of the form 

= y^ = u2QX^^a''\ (15) 

In light of Eq. (13), we take Eq = M/(37ri/io). Finally, we measure time in units of the 
viscous time-scale at the inner edge, t — toT, to — rfn/j^io, and define 770 = 1^20/^^10- With 



^This expression for S/;at applies to a nonmagnetic disk. When magnetic fields thread the disk, the 
functional form of J{r) is modified in a model-dependent way (see Lai 1999, appendix A for examples). 
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these assumptions, Eqs. (7)-(9) simplify to 



da 1 d 
dr X dx 



1" 



0, (16) 



' 2 



ax{x'^fiy 



where a prime denotes d/dx. The dimensionless torque is given by 



n. 
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(18) 



where 



3 sin^ ^ 
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3.5 



5.4 



0.5 



-3.5 /sin^^ 



f. = |,--Ccos^^ = 2l(^) 



V 0.5 

3.5 



cos^ 6' 
0.5 



(19) 
(20) 



are dimensionless constants independent of r and t, parametrizing the problem. We will first 
set D{x) = 1 to make contact to the linear stability analysis in Shirakawa & Lai (2002b). 
Later in this paper, we will use Eq. (3) and set for concreteness if(r)/rj„ = 0.02 close to the 
inner edge, so that 

(21) 



D{x) = max (^Vx^ - 1 , 0.2^ . 



(The specific choice of D{x) infiuences our qualitative results only marginally, as the radial 
dependence of the precession torque in Eq. (1) is dominated by the piece.) 

While eqs. (16) and (17) can be applied for general viscosity law [eq. (15)], in our 
calculations we shall restrict to the cases where 1^2/^2 = ^^20/^^10 = Vo- The ratio rjo measures 
the disk's resistance to warping; we will consider different values of rjo, ranging from r]o = 1 
to ?7o > 1 ^• 



^For a Keplerian disk, assuming isotropic viscous stress parametrized by the a-ansatz, the ratio ^2/2^1 is 
of order l/(2a^) for linear warps (Papaloizou & Pringle 1983; Ogilvie 1999). The large 1^2/2^1 (for a ^ 1) is 
due to the fact that horizontal motions induced in the disk by the warp are resonantly driven. Such resonant 
behavior tends to be diminished in the nonlinear regime and also when deviation from Keplerian rotation is 
present (Ogilvie 1999). 
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Finally, we give the boundary conditions in dimensionless variables: 
di 

g^i^in) = 0, i{Xout) = hut (22) 

a{Xin) = (Tfiat{Xin), (j{Xout) = flat{Xout) , (23) 

where, by definition, Xin — 1, and afiat = ^fiat/^o follows from (13) and (15): 

aflat = [(1 - (1 - J^n)x-'^') X"^^] '^^'"""^ . (24) 

4. Numerical Results 

Equations (16) and (17) with boundary conditions given by Eqs. (22) and (23) are solved 
with a Crank- Nicholson scheme as described in the appendix. 

4.1. Evolution Into the Nonlinear Regime 

We first consider an almost fiat, unstable disk which we follow through the linear growth 
of the warping mode into the nonlinear regime. The particular parameters we use are 
flp = 10, — 10, J'in = 1, Tjo = 1, Pi = (^2 — 0.6, ai — a2 — and set D{x) = 1. The linear 
stability analysis for these parameters was presented in Shirakawa & Lai (2002b), where the 
disk was found to be unstable to the magnetic warping instability. As initial conditions for 
the evolution we take (7{x)\^^^ = (Jfiat{x), and perturb 1 by ~ 10~^ degrees away from a). 

Fig. 1 presents the tilt-angle at the inner edge, j3in = I3{xin) and the precession frequency 
at the inner edge, d'-fin/dr as a function of time r. The disk is indeed unstable and a growing 
mode emerges. At early times, the tilt-angle grows exponentially with e-folding time of 
^e-fold — {dlog Pin/ dr)'^ = 0.59 and precession frequency djin/dr = —4.22. These numbers 
agree very well with the hnear stabihty analysis of Shirakawa & Lai (2002b). At late times 
the growth slows down, and the tilt-angle saturates around P{rin) ~ 75°. Eventually, the disk 
settles into a solution where the whole disk precesses with uniform precession frequency, i.e. 
the tilt-angle profile /9(r) is independent of time, whereas the twist angle increases linearly 
in time, j{r, r) = 7o(?^) + ^^pT, with the precession frequency cUp ~ —0.72. 

Fig. 2 presents the tilt-angle P{r) as a function of radius at different times. One sees 
clearly that after the mode has saturated, the tilted region of the disk is larger than during 
the linear phase. The profiles of the tilt-angle provide an explanation for the reduction of 
precession frequency as the disk enters the non-linear regime: The precessional torque given 
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by the first term of Eq. (18), is suppressed by the factor I ■ u — cos/3. Furthermore, a 
more extended disk is involved, reducing the precession rate further. In the final steady- 
state, the precession period is P = 27i/\ujp\ = 8.7, while the warp diffusion time T^[q is of 
order \d\/dr\~^i'2^ ~ = x^''^/i'2o (in dimensionless units). The disk warp extents to 

^ ~ 3 (see Fig. 2) so that T^^iff ~ 4.7 < P. Thus the different parts of the warped disk 
can communicate with each other through viscous diffusion on a timescale shorter than the 
precession time, making rigid-body precession possible. Also, both P and Tjjff are much 
larger than the e-folding time T^^^^i^ 0.59 of the growth of linear perturbations. The 
inequalities Tg.foi(j < T^iQ < P hold for all evolutions presented in this paper. 

Fig. 3 shows the surface density relative to the flat disk, a /a fiat at different times. As the 
evolution proceeds, the disk gets depleted. The reason for this is that as the disk warps, the 
second term in Eq. (9) gives rise to increased advection relative to the flat disk. Although 
the disk-tilt is conflned to relatively small radii (e.g., /3 < 2° for x > 5; and /3 decreases 
exponentially with larger x), the surface density is depleted over a larger region with the 
deviation from (TfUit decreasing only as The outer disk depletes on the viscous time- 

scale of the outer disk; this explains the very slow relaxation to the steady-state in Fig. 1. 
The fact that a changes at large radii necessitates the use of a large computational domain, 
although the disk is only warped in the inner region. We typically place the outer boundary 
at Xout = 1000, which results in accuracy of about 1%. 



4.2. Nonlineeir Steady State of Warped Disk 

The growing mode is only a transient phenomenon — one will most likely observe only 
the nonlinear steady state solution. Therefore we now study the properties of the steady- 
state solution as a function of four of the parameters that characterize the system: f ^, flp, 
Jin and r/o (for deflniteness, the other parameters are flxed as in §4.1). 

First we consider the effect of varying F^^, (keeping Op = 10, J7i„ = 1 and ri^ = \ and 
D(,t) = 1 as in §4.1). For different F^, we perform evolutions until relaxation to steady- 
state, and then record the tilt angle (5in and the disk precession frequency 0;^ = d'^in/dr. 
Figure 4 presents the results of this computation. For f ^ > 6.2, the disk is unstable to the 
magnetic warping instability (in agreement with the linear stability analysis) ; with increasing 
torque F^, the steady-state tilt-angle increases. There is only a small window of torque 
parameters, for which the steady-state disk is moderately warped (say, < (5in < 60°). 
Hence, this computation suggests that, whenever the magnetic warping instability operates, 
it is likely that the inner disk will be signiflcantly tilted. 
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Figure 4 also shows the steady-state precession frequency Up as a function of Ty,. At the 
onset of warping instabihty (r^^ ~ 6.2), uj^ ^ —3.7. As increases, \u!p\ drops very rapidly: 
Doubling the warping torque decreases \ujp\ by roughly a factor of ten. The reason is that 
for larger F^^,, the disk tilt is larger (approaching 90° very closely) and the warping region is 
more extended (see Fig. 5). Note that the surface density of the F^ = 100 evolution deviates 
less from a fiat than the one for f ^ = 10; this is because the local disk-warp is smaller in 
the former case (91/9 In r < 0.24 as compared with 91/9 In r < 0.94), so that less additional 
advection is introduced by the warping. The strong dependence of ujp on F^ opens the 
possibility to account for a wide variety of low-frequency quasi-periodic oscillations with the 
nonlinear behavior of a warped disk driven by the magnetic warping instability (see §5). 

We now use eq. (21) for the function D{x) and extend the parameter searches to Qp, Jin 
and 770- Figure 6 presents results for different choices of these parameters. For each curve, 
evolutions were first performed for large until the steady state is reached. Then was 
gradually lowered so that a sequence of steady-state solutions was traced out. (When 
becomes smaller than a certain critical torque, the disks suddenly flattens completely; see 
§4.3.) 

Several general features of the steady-state disks appear to be robust: In all cases, Up 
strongly depends on F^. Furthermore, f3in becomes significant (> 60°) for Tyj even slightly 
exceeding the critical value for the onset of the warping instability. The precession frequency 
uOp depends strongly on Op as well, as comparison between the Op = 10 and Op = 25 curves 
reveals: For the same V^j, changing Op by a factor of 2.5 results in a change of uip by about 
one order of magnitude. 

Consider now the effect of larger viscosity ratio ryo = ^^20/^^10 [= 12.5 and 50, corre- 
sponding to a = 0.2 and 0.1 in ^2/^1 = l/(2ct^); see §3.2]. A larger ^2/^1 represents a disk 
more resistive to bending, so that a larger warping torque F^^, is required for the magnetic 
warping instability to operate. Also, (3in approaches 90° more slowly as f ^ is increased. 
Figure 7 shows the steady-state disk profiles for the case of 770 = 50. These differ appreciably 
from the profiles obtained with 770 = 1 (as shown in Fig. 5): The tilted region of the disk 
extends to larger radii {x ^ 20) for 770 = 50 than for 770 = 1; the disk is also locally fiatter 
(91/91nr < 0.13) and more depleted (with a/ufiat being as small as 0.03). 

The inner boundary condition is somewhat critical. The disk behavior described above 
works for relatively large Jin (^ 0.2; depending on the other parameters). However, reducing 
Jin further leads to singularities in the disk: For evolutions with smaller Jin, say, Jin — 0.1, 
it appears that the surface density tends to zero at some finite radius Xging away from the 
inner edge. The tilt angle (3 is non-zero inside Xsing and close to zero outside Xsing- Overall, 
it appears that the disk pinches off and two regions remain: a flat outer disk with inner 
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boundary at Xsing, and a strongly tilted, rapidly precessing inner disk inside Xsing (which 
may be viewed as part of the magnetosphere) . 

4.3. Discontinuous Behavior: Hystereses 

In Figure 6, the tilt angle along each curve does not approach zero continuously 
reduced, but exhibits a discontinuity. A closer examination reveals a hysteretic 
behavior of warp disks. 

Consider the 770 = 12.5 case in Fig. 6 as an example (see Fig. 8). If we start with a slightly 
perturbed disk {(3 <^ 1, a — (Jfiat), then the perturbation decays for < Tcriti ~ 31.3 and 
the disk settles down into the flat state. For Vyj shghtly larger than Tcriti: the growth rate of 
the perturbation is proportional to F^ — Fcriti; nonetheless the perturbation grows to large 

amplitude (with /9.j„ reaching 70° in the final steady state). Thus, starting from small 
perturbations, the disk evolves into a steady state which depends discontinuously on (see 
the solid line of Fig. 8). 

The situation is different if one starts from a highly warped state: We first evolve the 
disk with > Tcrui until the final (strongly warped) steady state is reached. We then 
reduce Ty^ gradually; for each new value of F^ we evolve the disk to a new steady state 
(using the "previous" steady-state solution as the initial condition). As long as F^ > Tcrit,i 
this procedure produces the same sequence of steady-state disks we found when starting 
from a slightly perturbed disk. However, we find that even below Tcrit,i, there still exits 
steady-state warped disks (see the dashed curve in Fig. 8). We can follow this sequence 
down to a second critical value Tcrit,2 ~ 22.5. Slightly above Tcrit,2, the tilt angle depends 

— — — 1 /9 

on as Pin — Pcru + 0, {Vw — ^crit,2) for somc constant a and (5crit ~ 50°. As we reduce 
Ty, below Tcrit,2-i this warped steady-state sequence terminates abruptly and the disk relaxes 
to the fiat state. 

Thus, for rcrit,2 < < Fcrit,!, both a fiat disk and a warped steady-state disk solution 
exist. The steady state for such a disk depends on the history of the disk. When the disk 
parameters vary slightly — for example because the accretion rate varies — the disk might 
switch discontinuously between the two states. As the disk moves across a critical value, 
a strongly warped disk might suddenly "turn off" and become flat, or a flat disk might 
suddenly become unstable with the warp growing to signiflcant values. Since the local disk 
warp is fairly small at each radius (the maximum |9i/(?lnr| being approximately 0.3), it 
seems unlikely that the hysteresis is merely an artifact of the adopted simplifying evolution 
equations (see §3). 
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We find that all parameter sets examined in Fig. 6 exhibit such a hysteresis (e.g., 

for flp = 10, Jin = Tjo = 1: Tcrit,2 ~ 12.5,rcrit,i ~ 16.9; for = 25, Jin = r]o = 1: 
^crit,2 ~ 18,Tcrit,i ~ 38.2; for Op = 10, Jin = 0.5, i]o = 1: rcrit,2 ~ 15,rarit,i ~ 20.3; for 
Op = 10, Jin = l,Vo = 50: r^it,2 ~ 42.5, T^rit,i ~ 69.7). 

5. Discussion and Conclusion 

In this paper we have demonstrated through numerical simulations that accretion disks 
around magnetic stars can develop into a warped steady-state in which the disk attains 
a fixed (warped) shape and precesses rigidly. This work extends our previous (local and 
global) linear analysis of magnetically driven warping and precession of accretion disks. We 
find that whenever the magnetic warping instability criterion is satisfied in the linear regime, 
the inner disk (close to the magnctosphere boundary) is likely to be highly warped. The 
steady-state precession frequency spans a wide range, depending sensitively on the warping 
and precession torques (and hence on the parameters of the system). 

5.1. Limitations of the Models 

Before drawing any astrophysical conclusion of our work, it is useful to recall some of the 
limitations of our models. First, the nonlinear warp equations we adopted in our simulations 
are based on phenomenological descriptions of viscosities (see §3). Second, the magnetic 
torques formulae were derived under several assumptions/ansatz: e.g., locally flat disk, the 
use of the free parameter C to parametrize the magnetic field twist in quasi-steady state (see 
§3). Finally, there are intrinsic uncertainties associated with the inner boundary conditions 
of the disk. In our problem, the inner boundary is located at the magnctosphere radius, 
where the transition between a Keplerian disk and a corotating magnctosphere occurs. The 
physics that determines this transition is obviously complicated (see references cited in §1). 
The magnetic warping and precession torques are steep functions of radius and are maximal 
close to the inner disk edge. Therefore the details of our numerical results depend sensitively 
on the physics at the inner radius of the accretion disk. In our calculations we adopted the 
simplest inner boundary conditions for the disk. A more rigorous resolution of this problem 
will have to involve studying the coupling of the warped disk and the magnctosphere. 

Given these uncertainties and limitations, we cannot be sure whether some of the fea- 
tures we found for the warped disks (such the hysteresis behavior discussed in §4.3) corre- 
spond to reality in any way or simply represent artifacts of our models. Nevertheless, our 
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numerical experiments show that the general behaviors of magnetically driven warped disks 
that we found (as summarized at the beginning of §5) are robust. 



5.2. Astrophysical Implications 

The most important feature of magnetically driven disk warping studied in this paper 
is that the disk naturally evolves into a steady state which is highly warped near the inner 
radius (the magnetosphere boundary). This is in contrast to other warping mechanisms (e.g., 
those due to radiation pressure or irradiation-driven disk wind) which operate from outside- 
in. Such a highly warped disk can lead to modulation of the observed radiation flux either by 
obscuring the central star or by changing the reprocessed disk emission. As mentioned in §1, 
there is growing observational evidence for quasi-periodic oscillations (QPOs) or variabilities 
in radiation fluxes induced by warped inner disks in various systems, ranging from YSOs 
(e.g. AA Tauri; see Bouvier et al. 1999) to X-ray binaries (e.g. Her X-1: see Deeter et 
al. 1998; Scott et al. 2000; see also Shirakawa & Lai 2002b for a hst of other X-ray binaries 
showing mHz QPOs). 

Another interesting feature of our models is that the disk can be either in a completely 
flat state (warping stable) or in a highly warped state. Thus for a given system, as the 
disk parameters (e.g. accretion rate) vary, the disk can switch abruptly between the two 
states, leading to the appearance/disappearance of QPOs (with the corresponding change 
in radiation flux or spectrum). Indeed, in X-ray binaries, there are many examples where 
QPOs occur only in certain spectral states and not in the others (e.g., van der Khs 2000). 

Our calculations showed that the steady-state precession frequency of the disk can be 
much smaller (by up to several orders of magnitude, depending on the system parameters, 
cf. Fig. 6) than the frequency of the linear mode; the latter was approximately equal to 
the precession frequency at the disk inner radius (see Shirakawa & Lai 2002a,b for typical 
numbers for X-ray binaries). Thus, for magnetized neutron stars (with surface magnetic field 
~ 10^^ G), the global precession frequency is 

j^^^^Mj^ = -(1.2 mRz)Ai^l^M-l/\;'^/^J:j'D-^sin''e, (25) 

where the dimensionless parameter A can be much smaller than unity, and = ///(10^° G cm^), 
M1.4 = M/(1.4M0), E4 = E(r,,)/(10^gcm-2), withrg = r,J{10' cm) = SArj 1^%' M^'/ M-^^' 
as well as M17 = M/(10^^gs~^). For typical parameters appropriate for T Tauri stars, the 
precession period is given by 
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where B — ji/R^ — lO^Ss G is the surface magnetic field of the star, R — {2Rq)R2 is 
the stellar radius, and Mi = M/(l Mq). Therefore magnetically driven precession can 
potentially explain some long-period QPOs/variabilities in X-ray binaries and YSOs^. It is 
of interest to note that the magnetically driven precession is retrograde, as observed in some 
systems (e.g.. Her X-1). 

Finally, if the outflows from YSOs arc produced from interaction between the stellar 
magnetic field and the disk, as in the X-wind models (Shu et al. 1994), then the magnetic 
effects studied in this paper may be responsible for the jet precession observed in several 
systems [sec Tcrquem et al. 1999 and references therein; see also Lai (2003) for possible 
warping instability in magnetically driven disk outflows]. 



This work has been supported in part by NSF Grants AST-9986740, AST-0307252, 
PHY-9900672 and NASA grant NAG 5-12034. 



A. Numerical method 



Equations (16) and (17) are essentially of the form 

da ^"^a — (J CAD 

dr dx^ dx 
_ jj9H _ ^dl _ ^ 

dr dx"^ dx 

with coefficients A,...,F depending nonhnearly on the variables a and 1 (for ease of notation, 
we omit the hat on 1 here). We discretize these equations using the Crank-Nicholson method 
(Press et al. 1992). Time-derivatives are discretized as usual, for example 

The index i labels the spatial grid-points, unbarred quantities like CTj denote values at the 
current time Tq (which are known), and barred quantities denote values at the new time 



^We have only included the magnetic torques in our calculations. In real astrophysical systems, there 
may be other torques that contribute to (or even dominate) the disk precession; for example the tidal torque 
in binary systems (e.g. Terquem ct al. 1999; Ogilvie & Dubus 2001) and the torque associated with general 
relativistic frame dragging effect in low-mass X-ray binaries (see Shirakawa & Lai 2002a). Note that even in 
those systems where other torques dominate disk precession, the magnetic warping torque studied here can 
still be important in exciting/maintaining the disk tilt. 
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To + At which are to be determined. Spatial derivatives are averaged over tq and tq + At, 
for example 

dxy.^ 2 (Ax)2 2 {AxY ■ ^ ^ 

Here, Ai is computed using the known values at tq, whereas Ai is a function of the desired 
values at Tq + Ar. This discretization leads to a nonlinear system of equations, which is 
solved iteratively: We start with a guess for a and 1 at time Tq + Ar, which we denote with 
tildes, 5"i, li- Based on these guesses, we compute coefficients Aj, . . . Fj. Now, we discretize 
spatial derivatives hke (A4) using Ai instead of Af. 

Aj q-j+i - 2ai + (7j_i Aj aj+i - 2ai + aj^i 
" 2 (Ax)2 2 (Ax)2 ■ ^ 

Substituting these difference expressions into Eqs. (Al) and (A2) now yields /mear tridiagonal 
systems of equations for a, and 1,. Solving these yields improved values {ai,\i) for the 
variables at the next time step which are used in place of {ai,li). We iterate this process 
until convergence. 

This scheme is computationally more expensive per time-step than an explicit method. 
However, it is second order accurate in time and unconditionally stable, so that one can 
take very large time-steps. In practice, we often exceed Courant-factors of 1000. These large 
time-steps are especially important during the slow relaxation toward a coherently rotating 
steady-state disk. In such a case, one is primarily interested in the time-independent final 
state, whereas the transient evohition toward this state is less important. 

The code uses an adaptive time step: Ar is adjusted such that the difference in so- 
lution at Tq and Tq + Ar is smaller than some threshold, typically 10^"^. Furthermore, the 
code can be used with different radial distributions of grid points, for example linearly or 
logarithmically spaced grid points. This is done via a mapping x ^ x — x{x), which re- 
lates a "computational" coordinate x (in which the grid is uniform) to the "physical" radial 
coordinate x. Rewriting "physical" derivatives in terms of the computational coordinate x, 

d dx d f dx\^ d'^ d'^x d 



+ 



dx dx dx ' dx"^ \dx J dx'^ dx'^ dx ' 

shows that the coefficients A, . . . F will be modified by this mapping, however, the principal 
structure of Eqs. (Al) and (A2) does not change. Usually, logarithmically spaced grid-points 
are used which accurately resolve the fine structure close to the inner edge, while allowing 
to move the outer boundary far out, typically to Xout = 1000. 
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We verified that our code is second order convergent in space and time and recovers 
the standard flat disk solution, Eq. (13), in the absence of warping/precession torques. 
Section 4.1 confirms that in the linear regime our code recovers the linear stability analysis 
of Shirakawa & Lai (2002b). 
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100 1000 10000 

time T 

Fig. 1. — Tilt angle j3in at the inner edge and precession frequency at the inner edge during 
evolution into the nonlinear regime (parameters flp = Tyj = 10, J7in = ^, Vo = ^, Pi = = 
0.6, (Ji = (72 = 0.) The inserts show early times. 
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Fig. 2. — Tilt angle /3 as a function of radius at different times during the evolution presented 
in Fig. 1. The solid lines illustrate the growth of the warp; at r = 5400, the disk has 
essentially reached its final steady state. The dashed hne represents the shape during the 
linear regime (r = 6, rescaled to coincide with the steady state-solution at the inner edge), 
which agrees well with the eigenvector of the linear stability analysis of Shirakawa & Lai 
(2002b). 
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Fig. 3. — Surface density relative to the unperturbed, flat disk at different times during the 
evolution presented in Fig. 1. The dashed line represents our initial condition, (t/ a fiat = 1- 
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Fig. 4. — The tilt-angle at the inner edge of the disk, (3in, and the disk precession frequency 
\ujp\ as a function of warping torque parameter r^j,. The other parameters are fixed to 
Clp = 10, Jin = l,r]o = 1, Pi = ^2 = 0.6, cTi = 0-2 = and D{x) = 1. 
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Fig. 5. — The steady-state disk profiles for = 10 (solid lines) and = 100 (dashed 
lines), with the other parameters the same as in Fig. 4. The top panel shows the tilt-angle 
/5 as a function of radius, the bottom panel the surface density. Note that for F^^, = 100, the 
disk is tilted out to larger radii, but the surface density deviates less from that of the flat 
disk. 
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Fig. 6. — Similar to Fig. 4, but for different sets of parameters. The short dashed hnes 
correspond to a run with the same parameters as in Fig. 4, Clp — 10, = 1, 770 = 1, but 
with D{x) given by Eq. (21). Each one of the remaining runs is obtained by changing the 
value of one of these parameters, as labeled. 
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Fig. 7. — The steady-state disk profiles for f ^ = 50 (solid lines) and f = 100, both with 
rjo = 50, D{x) from eq. (21) and the other parameters as in Fig. 4. The top panel shows the 
tilt-angle /3 as a function of radius, the bottom panel the surface density (cf. Fig. 5. 
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Fig. 8. — Hystereses of warped disks. Plotted is the steady-state tilt angle at the inner disk 
edge as a function of the torque parameter f^. See §4.3 for details. 



